%% pde_stokes
clc;
clear;

index = 1;

syms x y z nu
switch index
    case 1
        u1 = sin(x).^2.*(sin(2*y).*sin(z).^2-sin(y).^2.*sin(2*z));
        u2 = sin(y).^2.*(sin(x).^2.*sin(2*z)-sin(2*x).*sin(z).^2);
        u3 = sin(z).^2.*(sin(2*x).*sin(y).^2-sin(x).^2.*sin(2*y));
        p = cos(x).*sin(y).*cos(z);
end

u1dx = diff(u1,x);
u1dy = diff(u1,y);
u1dz = diff(u1,z);
u1dxdx = diff(u1dx,x);
u1dydy = diff(u1dy,y);
u1dzdz = diff(u1dz,z);
u2dx = diff(u2,x);
u2dy = diff(u2,y);
u2dz = diff(u2,z);
u2dxdx = diff(u2dx,x);
u2dydy = diff(u2dy,y);
u2dzdz = diff(u2dz,z);
u3dx = diff(u3,x);
u3dy = diff(u3,y);
u3dz = diff(u3,z);
u3dxdx = diff(u3dx,x);
u3dydy = diff(u3dy,y);
u3dzdz = diff(u3dz,z);
pdx = diff(p,x);
pdy = diff(p,y);
pdz = diff(p,z);

f1 = -nu*(u1dxdx+u1dydy+u1dzdz) + pdx;
f2 = -nu*(u2dxdx+u2dydy+u2dzdz) + pdy;
f3 = -nu*(u3dxdx+u3dydy+u3dzdz) + pdz;

disp("===u1导数===");
matlabFunction(simplify(u1dx), "Vars", {x,y,z})
matlabFunction(simplify(u1dy), "Vars", {x,y,z})
matlabFunction(simplify(u1dz), "Vars", {x,y,z})
disp("===u2导数===");
matlabFunction(simplify(u2dx), "Vars", {x,y,z})
matlabFunction(simplify(u2dy), "Vars", {x,y,z})
matlabFunction(simplify(u2dz), "Vars", {x,y,z})
disp("===u3导数===");
matlabFunction(simplify(u3dx), "Vars", {x,y,z})
matlabFunction(simplify(u3dy), "Vars", {x,y,z})
matlabFunction(simplify(u3dz), "Vars", {x,y,z})
disp("===p导数===");
matlabFunction(simplify(pdx), "Vars", {x,y,z})
matlabFunction(simplify(pdy), "Vars", {x,y,z})
matlabFunction(simplify(pdz), "Vars", {x,y,z})
disp("===定常stokes方程===");
matlabFunction(simplify(f1), "Vars", {x,y,z,nu})
matlabFunction(simplify(f2), "Vars", {x,y,z,nu})
matlabFunction(simplify(f3), "Vars", {x,y,z,nu})